#!usr/bin/perl -w
=head1                           Bacterial-Scan_workflow
### Usage:
#	-fqlist <s> : Required Input
#	column1: ID
#	column2: fastq1 file
#	column3: fastq2 file
#	
#	-outdir <s> : Output Directory Must be absolute path!!!
#	              A file called "Fq.Trim.xls" will be generated in outdir
#	                            A shell script "Fq.Trim.sh" will also be generated in outdir
##   ...
#step1:	perl Bacterial-Scan.pl  raw.fq.list 
=cut
my $java="java";
use strict;
use warnings;
use Getopt::Long;
use File::Basename;
use FindBin qw($Bin $Script);
my ($fqlist,$outdir,$kmer,$help);
GetOptions(
		"fqlist:s"=>\$fqlist,
		"outdir:s"=>\$outdir,
		"help"=>\$help,
	  );
$outdir||=`pwd`;chomp $outdir;
die `pod2text $0` if (@ARGV < 1 || $help);
my ($line,@inf);
my $trim="01.QC";		mkdir $trim unless (-d $trim);
my $assem="02.Assembly";	mkdir $assem unless (-d $assem);
my $modify="03.Modify";		mkdir $modify unless (-d $modify);
my $anno="04.Annotation";	mkdir $anno unless (-d $anno);
my $circle="05.circle";		mkdir $circle unless (-d $circle);
my $plas="06.plasmid";  	mkdir $plas unless(-d $plas);
my $shell="shell";		mkdir $shell unless (-d $shell);							
my $result="result";		mkdir $result unless (-d $result);

my $rtrim="$outdir/$result/01.RawDataQC";        mkdir $rtrim unless (-d $rtrim);
my $rassem="$outdir/$result/02.Assembly";        mkdir $rassem unless (-d $rassem);
my $rpred="$outdir/$result/03.GenePredict";      mkdir $rpred unless (-d $rpred);
my $rncRNA="$outdir/$result/04.rRNA-tRNA";       mkdir $rncRNA unless (-d $rncRNA);
my $ranno="$outdir/$result/05.Annotation";       mkdir $ranno unless (-d $ranno);
my $rcirc="$outdir/$result/06.circle";       	 mkdir $rcirc unless (-d $rcirc);

my $gobal_sample_id;

open IN, "$ARGV[0]" or die "can not open file: $fqlist\n";	open RESULT,">$outdir/shell/S02.final.sh";
open OUT,">$outdir/$trim/clean.fq.list";	open SCAF,">$outdir/$modify/scaf.summ.lst";	open GENE,">$outdir/$modify/gene.info.lst";
open NCR,">$outdir/$modify/ncRNA.lst"; open ANNO,">$outdir/anno.lst";
while($line=<IN>){
	chomp $line;	@inf=split /\s+/,$line;		next if($#inf!=2);	$gobal_sample_id=$inf[0];
	print SCAF "$inf[0]\t$outdir/$modify/$inf[0].scaf.fna.summary.xls\n";
	print GENE "$inf[0]\t$outdir/$modify/$inf[0].GeneInfo.txt\n";
	print NCR  "$inf[0]\t$outdir/$modify/$inf[0].scaf.fna\t$outdir/$modify/$inf[0].tRNA.gff\t$outdir/$modify/$inf[0].rRNA.gff\n";
	print ANNO "$inf[0]\t$outdir/$modify/$inf[0].faa\t$outdir/$anno/$inf[0]/NR/$inf[0].vs.nr.table.xls\t$outdir/$anno/$inf[0]/GO/GO.list\t$outdir/$anno/$inf[0]/COG/$inf[0].cog.anno.xls\t$outdir/$anno/$inf[0]/KEGG/pathway.txt\t$outdir/$anno/$inf[0]/Swiss/$inf[0].vs.swiss.table.xls\t$outdir/$anno/$inf[0]/CAZy/$inf[0].anno.xls\t$outdir/$anno/$inf[0]/CARD/$inf[0].card.anno.xls\n";

	open SH, ">$outdir/$shell/S01.1.$inf[0].sh";

	print OUT "$inf[0]\t$outdir/$trim/$inf[0].clip.1.fq.gz\t$outdir/$trim/$inf[0].clip.2.fq.gz\n";

	print SH "$java -jar $Bin/Trimmomatic-0.33/trimmomatic-0.33.jar PE -threads 8 -phred33 $inf[1] $inf[2] $outdir/$trim/$inf[0].clip.1.fq.gz $outdir/$trim/$inf[0].single.R1.fastq.gz $outdir/$trim/$inf[0].clip.2.fq.gz $outdir/$trim/$inf[0].single.R2.fastq.gz ILLUMINACLIP:$Bin/Trimmomatic-0.33/adapters/merge.fa:2:30:10 SLIDINGWINDOW:4:15 MINLEN:75\n";
	print SH "rm $outdir/$trim/$inf[0].single.R1.fastq.gz $outdir/$trim/$inf[0].single.R2.fastq.gz\n";
	print SH "$Bin/FastQC/fastqc -o $outdir/$trim/ $inf[1] $inf[2]\n\n";

	`mkdir -p $outdir/$assem/$inf[0]`;	
	print SH "cd $outdir/$assem/$inf[0]/\n";
	print SH "/lustre/sdb/zengl/lib/anaconda3/envs/py2.7/bin/python $Bin/software/SPAdes-3.9.0/bin/spades.py -t 24 -m 150 --careful --disable-gzip-output -1 $outdir/$trim/$inf[0].clip.1.fq.gz -2 $outdir/$trim/$inf[0].clip.2.fq.gz -o $outdir/$assem/$inf[0]/\n";
	print SH "mv scaffolds.fasta ../$inf[0].scaffolds.fasta\n";
	print SH "perl $Bin/Bacteria_Scan/filter_spa_scaff.pl ../$inf[0].scaffolds.fasta > ../$inf[0].scaffolds.fasta.fil\n";
	print SH "perl $Bin/deal_fa.pl ../$inf[0].scaffolds.fasta.fil -format 2 -scafid $inf[0]_scaffold > $inf[0].scaf.fna1\n";	
	print SH "perl $Bin/deal_fa.pl $inf[0].scaf.fna1 -format 6 -len 500 >$inf[0].scaf.fna2\n";
	print SH "/lustre/sdb/linait/anaconda2/bin/python -m jcvi.apps.vecscreen mask --db=UniVec_Core $inf[0].scaf.fna2\n";
	print SH "cat $inf[0].UniVec_Core.bed |/lustre/sdb/linait/anaconda2/bin/maskFastaFromBed -fi $inf[0].scaf.fna2 -bed stdin -fo $inf[0].scaf.fna2.mask\n";
	print SH "/usr/bin/perl /lustre/sdb/linait/pipeline_micro/src/Assembly/Assembly_V2.6/lib/comm_bin/WGS_uplod_Seq.pl $inf[0].scaf.fna2.mask > $outdir/$modify/$inf[0].scaf.fna 2>$outdir/$modify/$inf[0].agp\n";
	print SH "rm -r $outdir/$assem/$inf[0]/\n\n";
	#print SH "/lustre/sdb/taoye/miniconda3/bin/abyss-pe name=$inf[0]-k64 k=64 in='$outdir/$trim/$inf[0].clip.1.fq.gz $outdir/$trim/$inf[0].clip.2.fq.gz'\n";
	#print SH "perl $Bin/Bacteria_Scan/assemle-summary.pl $outdir/$assem/$inf[0]/$inf[0]-k64-scaffolds.fa >$outdir/$assem/$inf[0]/$inf[0]-k64-abyss.scafSeq.summary;\n";

	`mkdir -p $outdir/$modify/$inf[0]`;
	open OOO,">$outdir/$modify/$inf[0]/$inf[0].ass.lst";
	print SH "cd $outdir/$modify/$inf[0]\n";
	print SH "perl $Bin/deal_fa.pl $outdir/$modify/$inf[0].scaf.fna -format 14 > $outdir/$modify/$inf[0].scaf.fna.summary.xls\n";
	print SH "perl /lustre/sdb/linait/pipeline_micro/bin/MicroGenome_pipeline3.3.pl --step 36 --SpeType B --contract_NO test --ass_seq $outdir/$modify/$inf[0]/$inf[0].ass.lst --step3 13 --ncRNA tRNA-rRNAd --outdir $outdir/$modify/$inf[0]/ --verbose\n";
	print SH "cp $outdir/$modify/$inf[0]/Result/03.Genome_Component/*/*/* ./\n";
	print SH "mv $inf[0].cds ../$inf[0].ffn ; mv $inf[0].pep ../$inf[0].faa; mv $inf[0].gff ../;mv $inf[0].tRNA.gff ../\n";
	print SH "mv $inf[0].denovo.rRNA.fa ../$inf[0].rRNA.fa;mv $inf[0].denovo.rRNA.gff ../$inf[0].rRNA.gff\n";
	
	print OOO "$inf[0]\t$outdir/$modify/$inf[0].scaf.fna\n";
	close OOO;
	print SH "ln -s ../$inf[0].scaf.fna\n";
	print SH "$Bin/software/Soap/soap2.21release/2bwt-builder $outdir/$modify/$inf[0]/$inf[0].scaf.fna\n";
	print SH "$Bin/software/Soap/soap2.21release/soap -a $outdir/$trim/$inf[0].clip.1.fq.gz -b $outdir/$trim/$inf[0].clip.2.fq.gz -D $outdir/$modify/$inf[0]/$inf[0].scaf.fna.index -o $inf[0].soap.pair.pe -2 $inf[0].soap.pair.se -l 35 -m 200 -x 800\n";
	print SH "$Bin/software/Soap/soap.coverage2.27 -cvg -il $outdir/$modify/$inf[0]/$inf[0].soap.list -refsingle $outdir/$modify/$inf[0]/$inf[0].scaf.fna -o $inf[0].coverage.txt -depthsingle ../$inf[0].coverage.fa\n";
	print SH "cd $outdir/$modify/\n";
	print SH "perl $Bin/software/GC_depth/gc_depth_20110530.pl $inf[0].scaf.fna $inf[0].coverage.fa -outdir $inf[0].GC_depth -windl 500 -y_mun 0,200,4 -x_mun 0,20,5\n";
	print SH "rm -r $outdir/$modify/$inf[0]/\n";
	print SH "perl $Bin/Bacteria_Scan/gene-statistic.pl $inf[0].scaf.fna  $inf[0].ffn > $inf[0].GeneInfo.txt\n\n";
	#print SH "$Bin/software/kmerfreq_v5.0/kmerfreq -k 17 -l $inf[0].fq2.list -L 300 >kmer.17.txt 2>kmer.17.log.txt\n";
	`mkdir -p $outdir/$anno/$inf[0]`;
	print SH "cd $outdir/$anno/$inf[0]\n";
	print SH "perl $Bin/produce_anno_sh.pl -name $inf[0] -db Bacteria $outdir/$modify/$inf[0].faa\n";
	print SH "sh Anno.sh\n\n";	

	print SH "cd $outdir/$circle/\n";
	print SH "perl $Bin/circle_figure/scaf_circle/plot_scaf_circle.pl $outdir/$anno/$inf[0]/COG/$inf[0].cog.anno.xls $outdir/$modify/$inf[0].gff $outdir/$modify/$inf[0].rRNA.gff $outdir/$modify/$inf[0].tRNA.gff $outdir/$modify/$inf[0].scaf.fna $inf[0].circle.pdf\n";
	open OUT1,">$outdir/$modify/$inf[0]/$inf[0].soap.list";	
	open CG, ">$outdir/$assem/$inf[0]/$inf[0].soap.config";
	print OUT1 "$inf[0].soap.pair.pe\n";
	
	print CG "max_rd_len=301\n[LIB]\n";
	print CG "avg_ins=450\nreverse_seq=0\nasm_flags=3\nrank=1\n";
	print CG "q1=$outdir/$trim/$inf[0].clip.1.fq.gz\nq2=$outdir/$trim/$inf[0].clip.2.fq.gz\n";
	#open PLAS,">$outdir/shell/S03.$inf[0].sh";
	#print PLAS "source activate plasflow\n";
	#print PLAS "cd $outdir/$plas\n";
	#print PLAS "python /lustre/sdb/taoye/miniconda3/envs/plasflow/bin/PlasFlow.py --input $outdir/$modify/$inf[0]/$inf[0].scaf.fna --output $inf[0].scaf.fna.p --threshold 0.9\n";
	#print PLAS "perl $Bin/filter_plasmid_result.pl $inf[0].scaf.fna.p > $inf[0].plasmid.summary.xls\n";
	#print PLAS "cat $inf[0].scaf.fna.p_plasmids.fasta > $inf[0].plasmids.fasta\n";
	#print PLAS "rm $inf[0].split*fa*\n";
	close SH;	close OUT1;    	close CG;	
	print RESULT "cp $outdir/$modify/$inf[0].GC_depth/gc_depth.png $rassem/$inf[0].gc_depth.png\n";
}

open SH3, ">$outdir/$shell/S01.2.stat.sh";
print SH3 "$java -jar $Bin/FastqStat.jar -i $ARGV[0]     > $outdir/$trim/rawdata.stat.main.xls\n";
print SH3 "$java -jar $Bin/FastqStat.jar -i $outdir/$trim/clean.fq.list > $outdir/$trim/cleandata.stat.main.xls\n";
print SH3 "cd $outdir/$trim/;ls *zip |while read dd;do unzip -x \$dd ;done;\n";
print SH3 "cp $gobal_sample_id*1.fq_fastqc/Images/per_base_quality.png          R1_per_base_quality.png\n";
print SH3 "cp $gobal_sample_id*1.fq_fastqc/Images/per_base_sequence_content.png R1_per_base_sequence_content.png\n";
print SH3 "cp $gobal_sample_id*2.fq_fastqc/Images/per_base_quality.png          R2_per_base_quality.png\n";
print SH3 "cp $gobal_sample_id*2.fq_fastqc/Images/per_base_sequence_content.png R2_per_base_sequence_content.png\n";

print SH3 "rm *html *zip\n";
print SH3 "perl $Bin/public/raw_clean2stat.pl $outdir/$trim/rawdata.stat.main.xls $outdir/$trim/cleandata.stat.main.xls >$outdir/$trim/data.stat.xls\n";
print SH3 "perl $Bin/Bacteria_Scan/st_scaf_summary.pl $outdir/$modify/scaf.summ.lst > $outdir/$modify/scaf.summary.xls\n";
print SH3 "perl $Bin/Bacteria_Scan/st_gene_info.pl    $outdir/$modify/gene.info.lst > $outdir/$modify/Gene.Info.xls\n";
print SH3 "perl $Bin/Bacteria_Scan/st_rRNA_tRNA.pl    $outdir/$modify/ncRNA.lst >$outdir/$modify/ncRNA.stat.xls\n";
print SH3 "perl $Bin/Anno/merge_anno_stat.pl $outdir/anno.lst >$outdir/$anno/annotation.stat.xls\n";
print SH3 "rm $outdir/$circle/*cmd.r\n";

close SH3;
open TABLE,">$outdir/table.stat.info.lst";
print TABLE "data_stat       $outdir/$trim/data.stat.xls\n";
print TABLE "assembly_stat   $outdir/$modify/scaf.summary.xls\n";
print TABLE "gene_stat       $outdir/$modify/Gene.Info.xls\n";
print TABLE "ncRNA_stat      $outdir/$modify/ncRNA.stat.xls\n";
print TABLE "annotation_stat $outdir/$anno/annotation.stat.xls\n";
close TABLE;
print RESULT "cd $rtrim/; ln -s $outdir/$trim/rawdata.stat.main.xls\n";
print RESULT "cd $rtrim/; ln -s $outdir/$trim/cleandata.stat.main.xls\n";
print RESULT "cp $outdir/$trim/*png $rtrim/\n";
print RESULT "cd $rassem/; ls $outdir/$modify/*fna     | while read dd;do ln -s \$dd ;done;\n";
print RESULT "cd $rpred/;  ls $outdir/$modify/*faa     | while read dd;do ln -s \$dd ;done;\n";
print RESULT "cd $rpred/;  ls $outdir/$modify/*ffn     | while read dd;do ln -s \$dd ;done;\n";
print RESULT "cd $rpred/;  ls $outdir/$modify/*gff     |grep -v \\rRNA  |grep -v \\tRNA | while read dd;do ln -s \$dd ;done;\n";
print RESULT "cd $rncRNA/; ls $outdir/$modify/*rRNA*   | while read dd;do ln -s \$dd ;done;\n";
print RESULT "cd $rncRNA/; ls $outdir/$modify/*tRNA*   | while read dd;do ln -s \$dd ;done;\n\n";
print RESULT "cd $ranno/;  cp -r $outdir/$anno/* ./\n";
print RESULT "cd $rcirc/;  cp $outdir/05.circle/*pdf ./\n";

print RESULT "cd $outdir/result/ ; cp -r $Bin/html_report/bact_scan/* .\n";
print RESULT "perl $Bin/html_report/table_html_height.pl $outdir/table.stat.info.lst >css/table.css\n";
print RESULT "cp $outdir/$trim/R1_per_base_quality.png          res/01.RawDataQC/per_base_quality.png\n";
print RESULT "cp $outdir/$trim/R1_per_base_sequence_content.png res/01.RawDataQC/per_base_sequence_content.png\n";

print RESULT "perl $Bin/html_report/table2html.pl $outdir/$trim/data.stat.xls >res/01.RawDataQC/data.stat.html\n";
print RESULT "perl $Bin/html_report/table2html.pl $outdir/$modify/scaf.summary.xls >res/02.Assembly/assembly.stat.html\n";
print RESULT "perl $Bin/html_report/table2html.pl $outdir/$modify/Gene.Info.xls    >res/03.GenePredict/gene.stat.html\n";
print RESULT "perl $Bin/html_report/table2html.pl $outdir/$modify/ncRNA.stat.xls   >res/04.rRNA-tRNA/ncRNA.stat.html\n";
print RESULT "perl $Bin/html_report/table2html.pl $outdir/$anno/annotation.stat.xls >res/05.Annotation/annotation.stat.html\n";
print RESULT "perl $Bin/html_report/anno2html.pl $outdir/$anno/$gobal_sample_id/NR/$gobal_sample_id.vs.nr.table.xls >res/05.Annotation/nr.anno.html\n";
print RESULT "perl $Bin/html_report/anno2html.pl $outdir/$anno/$gobal_sample_id/KEGG/pathway.txt >res/05.Annotation/kegg.anno.html\n";
print RESULT "perl $Bin/html_report/anno2html.pl $outdir/$anno/$gobal_sample_id/COG/$gobal_sample_id.cog.anno.xls >res/05.Annotation/COG.anno.html\n";
print RESULT "perl $Bin/html_report/anno2html.pl $outdir/$anno/$gobal_sample_id/Swiss/$gobal_sample_id.vs.swiss.table.xls >res/05.Annotation/swiss.anno.html\n";
print RESULT "perl $Bin/html_report/anno2html.pl $outdir/$anno/$gobal_sample_id/CAZy/$gobal_sample_id.anno.xls >res/05.Annotation/cazy.anno.html\n";
print RESULT "perl $Bin/html_report/anno2html.pl $outdir/$anno/$gobal_sample_id/CARD/$gobal_sample_id.card.anno.xls >res/05.Annotation/card.anno.html\n";

print RESULT "cp 02.Assembly/$gobal_sample_id.gc_depth.png                         res/02.Assembly/gc_depth.png\n";
print RESULT "ls $outdir/$anno/$gobal_sample_id/KEGG/pathways/*png |head -2 |tail -1 |while read d;do cp \$d res/05.Annotation/kegg.png; done;\n"; 
print RESULT "convert -density 150 -trim 05.Annotation/$gobal_sample_id/COG/COG.bar.pdf  res/05.Annotation/cog.png\n";
print RESULT "convert -density 150 -trim 05.Annotation/$gobal_sample_id/GO/GO.level2.pdf res/05.Annotation/go.png\n";
print RESULT "convert -density 150 -trim 06.circle/$gobal_sample_id.circle.pdf           res/06.circle/circle.png\n";

close IN; close OUT;close RESULT;
